/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 *
 * Redistribution and use in source and binary forms, with or without modification,
 * are permitted provided that the following conditions are met:
 * 
 * * Redistributions of source code must retain the above copyright notice, this 
 *   list of conditions and the following disclaimer.
 * * Redistributions in binary form must reproduce the above copyright notice, 
 *   this list of conditions and the following disclaimer in the documentation
 *   and/or other materials provided with the distribution.
 * * Neither the name of the dnAnalytics Project nor the names of its contributors
 *   may be used to endorse or promote products derived from this software without
 *   specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE 
 * DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
 * FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
 * DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
 * SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
 * CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
 * TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
 * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */

using System;
using System.Diagnostics;
using NUnit.Framework;
using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Solvers.Preconditioners;
using dnAnalytics.LinearAlgebra.Solvers;
using dnAnalytics.Math;
using System.Reflection;

namespace dnAnalytics.Tests.LinearAlgebra.Solvers.Preconditioners
{
    [TestFixture]
    [Category("Managed")]
    public sealed class IlutpPreconditionerTest : PreconditionerTest
    {
        private double mDropTolerance = 0.1;
        private double mFillLevel = 1.0;
        private double mPivotTolerance = 1.0;

        [SetUp]
        public void Setup()
        {
            mDropTolerance = 0.1;
            mFillLevel = 1.0;
            mPivotTolerance = 1.0;
        }

        private static T GetMethod<T>(Ilutp ilutp, string methodName)
        {
            var type = ilutp.GetType();
            var methodInfo = type.GetMethod(methodName,
                                            BindingFlags.Public | BindingFlags.NonPublic | BindingFlags.Instance | BindingFlags.Static,
                                            null,
                                            CallingConventions.Standard,
                                            new Type[0],
                                            null);
            var obj = methodInfo.Invoke(ilutp, null);
            return (T)obj;
        }

        private static SparseMatrix GetUpperTriangle(Ilutp ilutp)
        {
            return GetMethod<SparseMatrix>(ilutp, "UpperTriangle");
        }

        private static SparseMatrix GetLowerTriangle(Ilutp ilutp)
        {
            return GetMethod<SparseMatrix>(ilutp, "LowerTriangle");
        }

        private static int[] GetPivots(Ilutp ilutp)
        {
            return GetMethod<int[]>(ilutp, "Pivots");
        }

        private SparseMatrix CreateReverseUnitMatrix(int size)
        {
            SparseMatrix matrix = new SparseMatrix(size);
            for (int i = 0; i < size; i++)
            {
                matrix[i, size - 1 - i] = 2;
            }
            return matrix;
        }

        private Ilutp InternalCreatePreconditioner()
        {
            Ilutp result = new Ilutp();
            result.DropTolerance = mDropTolerance;
            result.FillLevel = mFillLevel;
            result.PivotTolerance = mPivotTolerance;
            return result;
        }

        internal override IPreConditioner CreatePreconditioner()
        {
            mPivotTolerance = 0;
            mDropTolerance = 0.0;
            mFillLevel = 100;
            return InternalCreatePreconditioner();
        }

        protected override void CheckResult(IPreConditioner preconditioner, SparseMatrix matrix, Vector vector, Vector result)
        {
            Assert.AreEqual(typeof(Ilutp), preconditioner.GetType(), "#01");
            Ilutp factorization = preconditioner as Ilutp;
            // Compute M * result = product
            // compare vector and product. Should be equal
            Vector product = new DenseVector(result.Count);
            matrix.Multiply(result, product);
            for (int i = 0; i < product.Count; i++)
            {
                Assert.IsTrue(Precision.EqualsWithinDecimalPlaces(vector[i],
                                                                  product[i],
                                                                  -Precision.Magnitude(s_Epsilon)),
                              "#02-" + i.ToString());
            }
        }

        [Test]
        public void SolveReturningOldVectorWithoutPivoting()
        {
            int size = 10;
            SparseMatrix newMatrix = CreateUnitMatrix(size);
            Vector vector = CreateStandardBcVector(size);
            // set the pivot tolerance to zero so we don't pivot
            mPivotTolerance = 0.0;
            mDropTolerance = 0.0;
            mFillLevel = 100;
            IPreConditioner preconditioner = CreatePreconditioner();
            preconditioner.Initialize(newMatrix);
            Vector result = new DenseVector(vector.Count);
            preconditioner.Approximate(vector, result);
            CheckResult(preconditioner, newMatrix, vector, result);
        }

        [Test]
        public void SolveReturningOldVectorWithPivoting()
        {
            int size = 10;
            SparseMatrix newMatrix = CreateUnitMatrix(size);
            Vector vector = CreateStandardBcVector(size);
            // Set the pivot tolerance to 1 so we always pivot (if necessary)
            mPivotTolerance = 1.0;
            mDropTolerance = 0.0;
            mFillLevel = 100;
            IPreConditioner preconditioner = CreatePreconditioner();
            preconditioner.Initialize(newMatrix);
            Vector result = new DenseVector(vector.Count);
            preconditioner.Approximate(vector, result);
            CheckResult(preconditioner, newMatrix, vector, result);
        }

        [Test]
        public void CompareWithOriginalDenseMatrixWithoutPivoting()
        {
            SparseMatrix sparseMatrix = new SparseMatrix(3);
            sparseMatrix[0, 0] = -1;
            sparseMatrix[0, 1] = 5;
            sparseMatrix[0, 2] = 6;
            sparseMatrix[1, 0] = 3;
            sparseMatrix[1, 1] = -6;
            sparseMatrix[1, 2] = 1;
            sparseMatrix[2, 0] = 6;
            sparseMatrix[2, 1] = 8;
            sparseMatrix[2, 2] = 9;
            Ilutp ilu = new Ilutp();
            ilu.PivotTolerance = 0.0;
            ilu.DropTolerance = 0;
            ilu.FillLevel = 10;
            ilu.Initialize(sparseMatrix);
            SparseMatrix l = GetLowerTriangle(ilu);
            // Assert l is lower triagonal
            for (int i = 0; i < l.Rows; i++)
            {
                for (int j = i + 1; j < l.Rows; j++)
                {
                    Assert.IsTrue(Precision.EqualsWithinDecimalPlaces(0.0,
                                                                      l[i,j],
                                                                      -Precision.Magnitude(s_Epsilon)),
                                  "#01-" + i.ToString() + "-" + j.ToString());
                }
            }
            SparseMatrix u = GetUpperTriangle(ilu);
            // Assert u is upper triagonal
            for (int i = 0; i < u.Rows; i++)
            {
                for (int j = 0; j < i; j++)
                {
                    Assert.IsTrue(Precision.EqualsWithinDecimalPlaces(0.0,
                                                                      u[i,j],
                                                                      -Precision.Magnitude(s_Epsilon)),
                                  "#02-" + i.ToString() + "-" + j.ToString());
                }
            }
            SparseMatrix original = l.Multiply(u) as SparseMatrix;
            for (int i = 0; i < sparseMatrix.Rows; i++)
            {
                for (int j = 0; j < sparseMatrix.Columns; j++)
                {
                    Assert.IsTrue(Precision.EqualsWithinDecimalPlaces(sparseMatrix[i,j],
                                                                      original[i, j],
                                                                      -Precision.Magnitude(s_Epsilon)),
                                  "#03-" + i.ToString() + "-" + j.ToString());
                }
            }
        }

        [Test]
        public void CompareWithOriginalDenseMatrixWithPivoting()
        {
            SparseMatrix sparseMatrix = new SparseMatrix(3);
            sparseMatrix[0, 0] = -1;
            sparseMatrix[0, 1] = 5;
            sparseMatrix[0, 2] = 6;
            sparseMatrix[1, 0] = 3;
            sparseMatrix[1, 1] = -6;
            sparseMatrix[1, 2] = 1;
            sparseMatrix[2, 0] = 6;
            sparseMatrix[2, 1] = 8;
            sparseMatrix[2, 2] = 9;
            Ilutp ilu = new Ilutp();
            ilu.PivotTolerance = 1.0;
            ilu.DropTolerance = 0;
            ilu.FillLevel = 10;
            ilu.Initialize(sparseMatrix);
            SparseMatrix l = GetLowerTriangle(ilu);
            SparseMatrix u = GetUpperTriangle(ilu);
            int[] pivots = GetPivots(ilu);
            SparseMatrix p = new SparseMatrix(l.Rows);
            for (int i = 0; i < p.Rows; i++)
            {
                p[i, pivots[i]] = 1.0;
            }
            SparseMatrix temp = l.Multiply(u) as SparseMatrix;
            SparseMatrix original = temp.Multiply(p) as SparseMatrix;
            for (int i = 0; i < sparseMatrix.Rows; i++)
            {
                for (int j = 0; j < sparseMatrix.Columns; j++)
                {
                    Assert.IsTrue(Precision.EqualsWithinDecimalPlaces(sparseMatrix[i, j],
                                                                      original[i, j],
                                                                      -Precision.Magnitude(s_Epsilon)),
                                  "#01-" + i.ToString() + "-" + j.ToString());
                }
            }
        }

        [Test]
        public void SolveWithPivoting()
        {
            int size = 10;
            SparseMatrix newMatrix = CreateReverseUnitMatrix(size);
            Vector vector = CreateStandardBcVector(size);
            Ilutp preconditioner = new Ilutp();
            preconditioner.PivotTolerance = 1.0;
            preconditioner.DropTolerance = 0;
            preconditioner.FillLevel = 10;
            preconditioner.Initialize(newMatrix);
            Vector result = new DenseVector(vector.Count);
            preconditioner.Approximate(vector, result);
            CheckResult(preconditioner, newMatrix, vector, result);
        }
    }
}
